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Drag Law of Two Dimensional Granular Fluids 

Satoshi Takada* and Hisao Hayakawa^ 


ABSTRACT 

The drag force law acting on a moving circular disk in a two-dimensional granular medium 
is analyzed based on the discrete element method (DEM). It is remarkable that the drag force 
on the moving disk in moderate dense and pure two-dimensional granular medium can be well 
reproduced by a perfect fluid with separation from the surface of the tracer. A yield force, being 
independent of the moving speed of the disk, appears if a dry friction between the granular disks 
and the bottom plate exists. The perfect fluidity is violated in this case. The yield force and the 
drag force diverge at the jamming point. 

Keywords: drag law, DEM, perfect fluidity, jamming transition. 


INTRODUCTION 

The drag law of a moving object in a medium is one of the most fundamental 
characterizations of rheology. The drag force on an object surrounded by a viscous 
fluid is proportional to the moving speed, whose dependence is known as Stokes law, 
if the moving speed is low. The drag is proportional to a fractional power of the moving 
speed, when the speed is faster. Then the power exponent reaches 2 in the high speed 
limit (ILamb 19451 IB atchelor 1967]) . 

The drag force on a moving object in a granular medium is completely different 
from that in the viscous fluid. Because granular materials behave as unusual solids 
and liquids (Jaeger et al. 19961, the drag force is thought to be a complex combina¬ 
tion of the force chains and fluid contribution. So far, some previous experiments 
reported that the drag force has non-trivial depth-dependence on the drag of cylinders 
(I Albert et al. 19991 IHill et al. 20051 IGuillard et al. 20131) and logarithmic dependence 
on the rotating frequency in a two-dimensional geometry and in a biaxial rotating cylin¬ 
der (Geng and Behringer 2004, Geng and Behringer 2005[ Reddy et al. 201 1| ). 

Most of these previous studies indicated the existence of two drag terms: one is 
independent of the moving speed, and another depends on that. Recently, Takehara and 
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her coworkers have performed a series of experiments to clarify the drag law acting on 
a circular disk in one layer granular medium by controlling the pulling speed V of the 
disk (ITakehara et al. 20T0t ITakehara and Qkumura 2014I) . They experimentally found 
that the drag force in a granular medium satisfies 


= Fo(c/>) + a(c/>)V^ 


( 1 ) 


where both a(<^) and Fo(<^) are proportional to (0c — with the area fraction 0 

and the jamming fraction 0c (ITakehara and Qkumura 20f4l) . It should be noted that Eq. 
([U) may be only valid for the moving object at relatively high speed as mentioned in 
their paper (ITakehara and Qkumura 20141) . though the expression itself can be extrap¬ 
olated ioV —j- 0. Some experiments on the impact of a hard projectile into granular 
beds also supported Eq. ([I]), though the density dependencies of the two terms are 
different ( [Katsuragi and Durian 2007 , Katsuragi and Durian 2013 Park et al. 20121 


Seguin et al. 2009). Although a previous study for three-dimensional simulation re¬ 


ported that the drag force is proportional to V (IHilton and Tordesillas 2013]) and the 

other studies (IChehata et al. 20031 [Wassgren et al. 2003^|Geng and Behringer 2O04||Geng and Behringer 2( 

Bharadwaj et al. 200^ suggested that the drag force logarithmically depends on V. 


Nevertheless, the drag force proportional to is natural for granular systems, which 


is consistent with Bagnold’s scaling (Bagnold 1954) and the direct impulse of granular 


particles hitting a moving object. The existence of the force being independent of the 
moving speed does not correspond to the drag law in a viscous fluid. To verify the va¬ 
lidity of Eq. ([U) and understand the mechanism to appear Eq. ([I]) are the central issues 
of this paper. 

Takehara and Okumura (2014) suggested that the existence of To(0) is related to 


the jamming transition ( 

Eiu and Nagel 1998, lO’Hem et al. 20021 lO’Hem et al. 20031 

Olsson and Teitel 2007llHatano 20081 Otsuki and Hayakawa 2009a 

Otsuki and Hayakawa 2009b 

Otsuki and Hayakawa 201 1[ Tighe et al. 2010 

ilNordstrom et al. 20101). It is. however. 


obvious that the jamming transition is unrelated to the existence of To(0)> because 
Fo(0) still exists far below the jamming density. The measurement of the drag force 
is relevant to know the viscosity in a viscous fluid. The viscosity of the granular fluid, 
however, obtained from the Couette flow is much larger than that obtained from the 
drag experiment. Indeed, an experiment for a granular jet ( [Cheng et al. 2007 ) as well as 
simulations (lEllowitz et al. 20131 IMiiller et al. 20141) suggested that the granular fluid 
can be approximately represented by a perfect fluid, though there exist counter ar¬ 
guments ([Sano and Hayakawa 2012[ [Sano and Hayakawa 2013]). Therefore, another 


purpose of this paper is to resolve the current confusing situation on the rheology of 
granular fluids. 


MODEL 

In this paper, we perform two-dimensional simulations in terms of the discrete 
element method (DEM) dCundall and Strack 19791) for a moving disk (the diameter D, 
the mass M, the position R, and the velocity V) surrounded by granular particles (the 
diameter di, the mass m*, the position r,, and the velocity Uj = r, for i-th grains and 
the number of grains is N) with or without the influence of dry friction characterized by 
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FIG. 1. (Color online) A schematic picture of our setup. We choose the 
puiiing direction of the tracer as negative x-direction. We aiso introduce 
the poiar coordinates (r, 9). Here, the arrow represents the moving direc¬ 
tion of the tracer. 


Coulombic friction constant /r between the bottom plate and the granular disks (Fig. 
[T])- Here, we consider both cases with and without the rotation of the disks and the 
tangential contacting forces. For the frictionless case, the equations of motion of the 
tracer which is the moving object in this study and the granular particle i are expressed 
as 

{mR= i^ex + -Fint(-R) ~ fiMgV , 

I miVi = Fint(rj) - firriigvi, 


where F]nt represents the interaction between grains satisfying 


= f 

Z—/ 7=0*' 1 


with/^ = Q{dij-rij){kn{dij-rij)rij-T]nrij-rij},andV = V/\V\ andvi = Vi/\vi\ 
are unit vectors parallel to V and vi, respectively. Here, Yl' denotes the summation 
under the condition j ^ i, and 0(x) is the step function, i.e. 0(x) = 1 for x > 0 and 
0(x) = 0 for otherwise. We characterize the tracer by i = 0, do = D, and ro = R. 
We also use dij = {di + dj)/2, the relative velocity Vij between i and j grains, and 

with the spring constant/c„ and the viscous parameter in the 
normal direction. It should be noted that Hertzian contact force in a two-dimensional 
system can be written as a term proportional to the compression with a logarithmic 
correction (I Johnson 19851 Gerl and Zippelius 1999[ Hayakawa and Kuninaka 2002l. 
In this paper, for simplicity, we adopt the linear spring model to represent the elastic 
force between contacting particles. 

For the frictional case, Eq. dH) still can be used with the replacement of the contact 
force by F-.yri) = Y'^=o(fPj + //j)’ where the tangential force fk is given by 


fij = ^{dij - y) min(/i,|/"|, ( 3 ) 

where we adopt /i^ = 0.2 for Coulombic friction constant between grains, fk = 

-kt^ij-T]tVt,ij, and tij = Here, Vt^ij = rij+rijx{diu:i+dju;j)-{rij-rij)rij 

with the angular velocity oji of z-th particle, ^ij is the tangential overlap vector between 
i-th and j-th particles defined by ^ij = dt'vt^ij{t') with the time to of first contact 
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of the particles, and kt and rjt are the spring constant and the viscous parameter in the 
tangential direction, respectively. We have introduced the function min(a, h) to select 
smaller one from a and b. In addition, an equation of motion for the rotation is given 
by 




(4) 


where d>i is the angular acceleration of i-th particle, J* = mid‘^/8 is the moment of 


inertia of i-th particle, Tj = ^{dij - x (/5 + flj)} is the torque of i-th 


particle, and Ctj is the rolling friction constant ( [Suyama et al. 2008] ). We assume that 
the mass density of each grain is identical and the system contains an equal number 
of two types of grains characterized by the diameters d and lAd to avoid the crys¬ 
tallization. It should be noted that the driving force i^ex = —Fex^x directly acts on 
the tracer i.e. moving disk, where is the unit vector in x-direction. The system 
reaches a steady state by the balance between T’ex and the other forces, Fint{R) and 
the dry friction force —jJiMgV between the tracer and the bottom plate. Thus, we 
obtain a steady motion of the moving disk to negative x-direction from the simula¬ 
tion of Eq. (I2l). We adopt the values of parameters rin = 0.75\/mfc„ corresponding 
to the restitution constant e = 0.92, = AA7md-^m/kn and the time increment 

At = 0.002^m/A:„, where m is the mass of the grain of the diameter d. The value 
of Cw is equal to Co; = 2Qmdy/d/g in terms of the gravitational acceleration g instead 
of kn- We have checked that the motion of the tracer is almost same in the range 
lOmdydd/g < (lj < 50md^/d/g when we perform the simulations using the identical 
initial condition. The values of kt and rft are chosen as kt = {2/7)kn and r]t = (2/7)?7„, 
respectively ([Thompson and Grest 1991|) and k^ is chosen as = SOOm^f/d with the 


gravitational acceleration g. For dry friction between the grains and the bottom plate, 
we examine /r = 0.001 and /r = 0. The reason why we adopt such a small value 
for /i is that it is difficult to obtain a steady motion for a wide range of the external 
force for larger /r. We adopt the velocity Verlet algorithm for the time integration of 
the equations of motion. The system size is basically fixed to be 210d x 105d. As a 
result, the number of grains N depends on the area fraction, though a typical value is 
A ~ 2 X 10^. The grains are located at random without any motion and overlap at time 
t = 0. For simplicity, we adopt the periodic boundary condition for x-direction and 
we assume that the boundary in ^/-direction is composed of particles whose mass and 
curvature are infinite with the spring constants kn, kt and the viscous parameters rjn, 
rjt, respectively. We have two control parameters, the area fraction 0 where the tracer 
is not included and the external force Fp,,. 


RESULTS 

System without dry friction 

First, we analyze a frictionless or a frictional system consisting of a large tracer 
disk and a collection of granular disks without the influence of the dry friction i.e. 
yU = 0 in Eq. (|2l). Figure [2] displays the density profile obtained from our DEM for 
Fex = 0.2knd and 0 = 0.76, and the streamlines of both our DEM (red solid lines) and 
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FIG. 2. (Color online) The density profiie (coior scaie) and the stream- 
iines (the red soiid iines) obtained from our DEM of frictionai grains for 
(j) = 0.76, Fex = 0.2knd, and the streamiines of the perfect fiuid (the biack 
dashed iines), where we have used D = lOd. Here the fiow direction is 
from ieft to right in the frame that the tracer is stationary. 



the perfect fluid (dashed lines) in the frame that the tracer is stationary. The streamlines 
are obtained by averaging over time during 100 y/m//c„ and 10 ensembles. It is notable 
that the density is almost uniform except for the cavity right behind the tracer, where 
no grains exist in this region. The streamlines obtained by the DEM are smooth and do 
not contain any vortex excitation. It is remarkable that the streamlines of the DEM are 
well reproduced by those of the perfect fluid except for that in the cavity. We stress that 
there is no contribution to the drag force from the cavity, because there are no grains 
colliding to the tracer there. 

Eigures[3];a) and (c) represent the velocity fields around the tracer against the polar 
angle 9 for Fex = 0.2knd and 0 = 0.76 for frictionless grains and frictional grains, 
respectively. These figures clearly support that the radial component of the granular 
flow on the surface of the tracer is almost zero, and the polar component of the granular 
flow can be approximately represented by a sinusoidal function of 6^ for 6*o < 6^ < 
27r — 9o and 0 for —Oq < 9 < 9q. Here, the separation angle, 9q (which is nearly 
equal to 80° for frictionless disks and 70° for frictional disks), is almost independent 
of V and 0 (Eigs. Ob) and (d)). If 6*o = 7r/2, it can be interpreted as the result of 
direct impulses of colliding grains. However, 9o is a little smaller, as seen in Eigs. Ob) 
and (d). The mechanism to have smaller 9o might be understood by the finite granular 
temperature effect. 

Eet us examine the drag force of the perfect fluid with the separation angle 6*o. It is 
well known that the pressure around a cylinder in the perfect fluid is given by 

p = Poc + ^V‘^(l-4cos^9), (5) 

for an irrotational incompressible perfect fluid, where p is the mass density of the 
granular fluid (IBatchelor 19671) . Because the far field pressure Poo is the impulse per 
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FIG. 3. (Color online) (a) The velocity fieid around the tracer against 9 
for Fex = 0.2knd and 0 = 0.76 for the frictioniess disks where we piot the 
radiai component Ur (the red open circies) and the poiar component ue 
(the biue open squares), (b) The veiocity dependence on the separation 
angie 6o for the frictioniess disks, where the separation angie is defined 
by the angie where ug deviates from a sinusoidai function, (c) The veiocity 
fieid against 6 for the frictionai disks, (d) The veiocity dependence on 6o 
for the frictionai disks. The dashed iine is the average over the resuits. 
Here, we have introduced dimensioniess quantities u* = Ur^/m/kn/d, Ug = 

ug\/m/kn/d and y* = V^mjkjd. 


unit cross line at the boundary, we adopt the expression p^o = {I + e)pV‘^, where e is 
the restitution eonstant. From the integration of the pressure aeting on the surfaee of 
the tracer, Fdrag = — Jg^~^°{D/‘2)d0p cos 9, we obtain the drag foree 


.^drag 


3 + 2e 


- sin^ 9o sin 9oDpV^. 


( 6 ) 


Note that F^rag should be zero if there is no separation, 6*o = 0, in Eq. db]). The granular 
fluid, however, has a finite eontribution even if the viseosity is zero beeause of the sepa¬ 
ration of the flow dSouthwell and Vaisey 19461). The expression dH) indieates that there 
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FIG. 4. (Color online) The dimensionless drag force -^drag — ^dvag/knd 
divided by 4> against V* for (a) the frictioniess disks and (b) the fric- 
tionai disks, respectiveiy, for three densities: </> = 0.72 (red open circies), 
0 = 0.74 (biue open squares), and 0 = 0.76 (pink open triangies), with 
Fq = Fo/knd. The dashed iines are theoreticai ones obtained from the 
perfect fiuid where the separation angies are, respectiveiy, determined 
by Figs.Htb) andHtd) for frictioniess and frictionai grains. 


is no yield force for pure two-dimensional cases, in contrast to previous experimental 
results (ITakehara et al. 20T0t ITakehara and Okumura 20 111) . It is astonishing that this 
simple formula dH) well reproduces the result of our simulation for 0 = 0.72,0.74 
and 0.76 without any fitting parameter (Fig. lU, while the formula dH) deviates from 
the simulation results when V* = V\/m/kn/d is larger than 0.3. This result is 
interesting because the perfect fluidity observed in granular jets or granular fluids 
(Cheng et al. 2007; lEllowitz et al. 20131 IMtiller et al. 20141: IBlumenfeld et al. 20101) is 
quantitatively verified in our setup, at least, for relatively slow flows and moderate 
dense granular medium. We should note that our problem can be converted into a jet 
problem for a circular target, if we use the frame of the stationary tracer. Because the 
flow has zero granular temperature at f = 0 , the excitation after the impact can be the 
origin of the viscosity. However, if the flow is slow and the target is circular, the ex¬ 
citation of the temperature by collisions is quite small and, thus, the flow in our setup 
can keep the perfect fluidity. We also note that Eq. db]) is no longer valid in the vicinity 
of the jamming point. 


The role of dry friction 

Now, let us consider the case of finite friction, i.e., /r 7 ^ 0 between the bottom plate 
and grains. Eor /r 7 ^ 0, the gravitational acceleration g can produce a new time scale 
\fdfg. Therefore, we expect that the yield force Fq is finite for /r 7 ^ 0. EigureOa) 
exhibits the result of our DEM for /r = 10“^. In this case, however, the perfect fluidity 
is violated. This violation might be related to the existence of force chains as in Eig. 
[5tb). Namely, the motion of grains is correlated with each other through the force 
chains. 
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FIG. 5. (Color online) (a) The relationship between and y* for fric- 
tioniess disks in 0 = 0.72 (red open circies), 0 = 0.74 (biue open squares) 
and 0 = 0.76 (pink open triangies). (b) A snapshot of the force chains 
(biue soiid iines) for frictioniess disks at 0 = 0.82, = 3.0, where iine 

width represents the strength of the interaction between contacting par- 
ticies. 


As the density increases, the drag force increases. Then it diverges at the jamming 
point. It seems that Eq. ([T]) is still valid, at least, for 0.1 < 17* < 0.3. In other words, 
we cannot fit the data by neither a logarithmic function nor a linear function of V for 
the drag. We determine the coefficient a in Eq. ([T]) by fitting the data in the range 
0.1 < 17* < 0.3, where Fq is estimated by extrapolating to E —)■ 0. Eigures [bta) and 
(b) show that a and Fq have the almost identical dependencies on the density as 

a(0) ~(0c - 0)"^, (V) 

i"o(0)~(0c-0)-^', (8) 


near the jamming point 0c where 0c = 0.8437, 0 = 0.277 ± 0.028, and 0' = 0.312 ± 
0.027, respectively (see Eigs.|^c) and (d)). It should be noted that the drag law deviates 
from the quadratic form for E* < 0.1. Similar behavior is also observed in a previ¬ 
ous experiment (lOkumura 20141) . We stress that the jamming point 0c in Eqs. (|7]) and 


which is a little larger than the previous estimations (Otsuki and Hayakawa 2009a: 


|Otsuki and Hayakawa 2012] ITakehara and Okumura 20141) . We also note that the re¬ 
sults cannot be represented by simple power laws as in Eqs. dV]) and ([8]) if we choose 
smaller 0c reported in the previous studies. The exponents of 0 and 0' in Eqs. (|7]) 
and dH]) are a little smaller than those by Takehara and Okumura (2014), though 0 is 
nearly equal to 0'. Therefore, we conclude that our simulation near the jamming point 
is qualitatively similar to that observed by Takehara and Okumura (2014) but quan¬ 
titative agreement is poor, which might be from the softness of the particles in our 
simulation. 


CONCLUSION 

In this paper, we performed two-dimensional OEMs to study the drag force acting 
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FIG. 6. (Color online) (a) The plot of divided by Fg against V* for fric- 
tioniess disks in 0 = 0.82 (red open circies), 0 = 0.83 (biue open squares), 
and 0 = 0.835 (pink open triangies). The dashed iine is a fitting curve 
in the range 0.1 < < 0.3 using Eq. ([1). (b) The piot of divided 

by (0c - 0)"^' against V* for 0 = 0.82, 0.83, and 0.835 with 0c = 0.8437 
and 0' = 0.312 ± 0.027. The dashed iine is a fitting curve in the range 
0.1 < 1^* < 0.3 using Eq. (H). (c) The density dependence of a, where 
dashed iine is given by q!( 0) ~ (0c-0)“^ with 0 = 0.277±0.028. (d) The den¬ 
sity dependence of Fq, where dashed iine is given by Fo(0) (0c -0) 


on the tracer for both frictionless and frictional granular disks with or without the 
influence of dry friction between the plate and grains. If there is no dry friction, we 
confirmed that the perfect fluid model with the separation of the flow can reproduce the 
quantitative behavior of the drag force for the moderate dense case such as 0 = 0.72, 
0.74, and 0.76. If there exists the dry friction, the yield stress Fo(0) appears and the 
perfect fluidity is no longer valid. In this case, the drag force and the yield force 
diverge at the jamming point, whose behavior is qualitatively similar to that observed 
in an experiment (iTakehara and Qkumura 20141) . but quantitative agreement between 
our results and their results is poor. 
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APPENDIX I. NOTATION 

The following symbols are used in this paper: 
d = diameter of smaller grains; 
do = diameter of the tracer, equivalent to 79; 
di = diameter of i-th particle; 
dij = average diameter of i-th and j-th particles; 

79 = diameter of the tracer; 

ffj = interaction between 7-th and j-th particles in the normal direction; 
f-j = interaction between 7-th and j-th particles in the tangential direction; 
fij = tangential interaction in the case of nonslip; 

Fq = yield force independent of the speed V; 

T^drag = force exerted on the tracer; 

Fex = external force exerted on the tracer; 

Fex = magnitude of F^x, 

Fint = interaction among the tracer and disks; 
g = gravitational acceleration; 
li = moment of inertia of 7-th particle; 
kn = spring constant in the normal direction; 
kt = spring constant in the tangential direction; 
nii = mass of 7-th particle; 

M = mass of the tracer; 

N = number of surrounding particles; 
r = radial coordinate in the polar coordinate; 

To = position vector of the tracer, equivalent to R, 

Vi = position vector of 7-th particle; 

Tij = distance between 7-th and j-th particles; 

Vij = unit vector parallel to relative position vector of 7-th and j-th particles; 
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iVi 


relative veloeity of i-\h and j-th partieles; 

position vector of the tracer; 

time; 

time of first contact of the particles; 
torque of i-th particle; 

unit vector parallel to the tangential force between i-th and j-th particles; 

radial component of the velocity near the tracer; 

radial component of the dimensionless velocity near the tracer; 

polar component of the velocity near the tracer; 

polar component of the dimensionless velocity near the tracer; 

unit vector parallel to the velocity of i-th particle; 

relative velocity between i-\h and j-th particles in the tangential direction; 

steady speed of the tracer; 

velocity of the tracer; 

dimensionless steady speed of the tracer; 

unit vector parallel to the velocity of the tracer; 

coefficient of the term proportional to the square of the speed V ; 

critical exponent of a with respect to the area fraction; 

critical exponent of the yield force Fq with respect to the area fraction; 

rolling friction constant; 

viscous parameter in the normal direction; 

viscous parameter in the tangential direction; 

angular coordinate in the polar coordinate; 

separation angle behind the tracer; 

The step function ©(x) = 1 and 0 for a; > 0 and a: < 0, respectively; 
Coulombic friction constant between the tracer and the bottom plate; 
Coulombic friction constant between grains; 
mass density of granular fluid; 

tangential overlap vector between i-th and j-th particles; 

area fraction of the system; 

jamming area fraction; 

angular velocity of f-th particle; and 

angular acceleration of f-th particle. 
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